Load libraries, load data, clean data

library(dplyr)
library(stringr)
library(lubridate)
library(Matrix)
library(igraph)
library(bipartite)
library(data.table)
library(ggplot2)
library(gganimate)
library(tidyr)
library(viridis)
library(asnipe)
set.seed(1234)

#Capture data from 2019/2020: refromat capture date as POISXct; retain PIT#, species, and capture date
pit_capture <-read.csv("All_PIT.csv", stringsAsFactors = FALSE, colClasses = "character")
pit_capture<- select(pit_capture, c(Lake, PIT, Length, Sex))
pit_capture$PIT<-as.character(pit_capture$PIT)
pit_capture<-pit_capture[!apply(pit_capture == "", 1, all),] #remove empty rows from initial read in
pit_capture<-pit_capture[which(pit_capture$Lake=="Parley"),]#limit to Lake Parley
head(pit_capture)
##      Lake    PIT Length Sex
## 2  Parley 497897    560   M
## 4  Parley 497857    720   F
## 6  Parley 497819    585   M
## 8  Parley 497887    690   M
## 10 Parley 497803    640   U
## 12 Parley 497862    740   U
length(unique(pit_capture$PIT))
## [1] 335
#Read in individual site data logs
s1 <-read.csv("SITE_1_2020.csv")
s2 <-read.csv("SITE_2_2020.csv")
s3 <-read.csv("SITE_3_2020.csv")
s4 <-read.csv("SITE_4_2020.csv")
s5 <-read.csv("SITE_5_2020.csv")
s6 <-read.csv("SITE_6_2020.csv")
s7 <-read.csv("SITE_7_2020.csv")
s8 <-read.csv("SITE_8_2020.csv")

all_sites <- rbind(s1,s2,s3,s4,s5,s6,s7,s8, Fill=NA)
# all_sites$PIT<-str_sub(all_sites$X, -6, -1) #note that in original data set length of some abbreviated pit numbers was only five characters (e.g #38305)
all_sites<-select(all_sites, c(Date, Time, PIT, Site))
all_sites$DATETIME<- as.POSIXct(paste(all_sites$Date,all_sites$Time), format='%m/%d/%y %I:%M:%S %p', tz="CDT") #Reformat as POSIXct
all_sites$PIT<-as.character(all_sites$PIT)
all_sites$Site<-as.character(all_sites$Site)
all_sites<-all_sites[-which(is.na(all_sites$PIT)),]
head(all_sites)
##      Date       Time    PIT Site            DATETIME
## 1 6/24/20 1:13:51 PM 567699    1 2020-06-24 13:13:51
## 2 6/24/20 1:13:52 PM 567699    1 2020-06-24 13:13:52
## 3 6/24/20 1:14:27 PM 567699    1 2020-06-24 13:14:27
## 4 6/24/20 1:14:32 PM 567699    1 2020-06-24 13:14:32
## 5 6/24/20 1:14:32 PM 567699    1 2020-06-24 13:14:32
## 6 6/24/20 1:14:34 PM 567699    1 2020-06-24 13:14:34
merged_all_sites<- merge(all_sites, pit_capture,  by="PIT")
head(merged_all_sites)
##     PIT   Date        Time Site            DATETIME   Lake Length Sex
## 1 38305 8/4/20 11:23:51 AM    4 2020-08-04 11:23:51 Parley    685   M
## 2 38305 8/4/20 11:23:50 AM    4 2020-08-04 11:23:50 Parley    685   M
## 3 38305 8/4/20 11:23:37 AM    4 2020-08-04 11:23:37 Parley    685   M
## 4 38305 8/4/20 11:23:49 AM    4 2020-08-04 11:23:49 Parley    685   M
## 5 38305 8/4/20 11:53:16 AM    4 2020-08-04 11:53:16 Parley    685   M
## 6 38305 8/4/20 11:53:06 AM    4 2020-08-04 11:53:06 Parley    685   M
merged_s1<-merged_all_sites[which(merged_all_sites$Site=="1"),]
merged_s2<-merged_all_sites[which(merged_all_sites$Site=="2"),]
merged_s3<-merged_all_sites[which(merged_all_sites$Site=="3"),]
merged_s4<-merged_all_sites[which(merged_all_sites$Site=="4"),]
merged_s5<-merged_all_sites[which(merged_all_sites$Site=="5"),]
merged_s6<-merged_all_sites[which(merged_all_sites$Site=="6"),]
merged_s7<-merged_all_sites[which(merged_all_sites$Site=="7"),]
merged_s8<-merged_all_sites[which(merged_all_sites$Site=="8"),]

#Site 1  
s1 <-read.csv("SITE_1_2020.csv") 
# s1$PIT<-str_sub(s1$X, -6, -1) #note that
# in original data set length of some abbreviated pit numbers was only five characters
s1<-select(s1, c(Date, Time, PIT, Site))
s1$DATETIME<- as.POSIXct(paste(s1$Date,s1$Time), format='%m/%d/%y %I:%M:%S %p', tz="CDT")
# #Reformat as POSIXct
s1$PIT<-as.character(s1$PIT) #toString(boatramp$PIT)
merged_s1_test<- merge(s1, pit_capture,  by="PIT")
length(unique(s1$PIT))
## [1] 88
length(unique(merged_s1$PIT))
## [1] 58

Data collection

In total, 170 unique individuals were recorded across 8 sites in Lake Parley the 2020 season.

 #how many uniquely tagged carp total and at each site?
(length(unique(merged_all_sites$PIT)))
## [1] 170
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='1')])))
## [1] 58
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='2')])))
## [1] 63
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='3')])))
## [1] 57
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='4')])))
## [1] 49
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='5')])))
## [1] 53
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='6')])))
## [1] 66
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='7')])))
## [1] 48
(length(unique(merged_all_sites$PIT[which(merged_all_sites$Site=='8')])))
## [1] 59
#' Define function to calculate number of visits per site
#' @param df- dataframe with $PIT ID column
#' @param num_days- number of days over which to loop
#' @param start_date- day at which to start looop
#' @return superfish df
calc_nvisit<-function(df,num_days,start_date){
  superfish<-data.frame(fishID=unique(df$PIT), nvisit=rep(NA,length(unique(df$PIT)))) 

for (i in 1:length(unique(df$PIT))){
    fishID<-unique(df$PIT)[i]
    visit_count<-0
    for (j in 1:num_days){
        sub_combined<-df %>% filter(DATETIME >= as.POSIXct(start_date+ (j-1)*86400, tz="CDT") & DATETIME <= as.POSIXct(start_date + j*86400, tz="CDT")) #86400 seconds/24 hours
        fishID %in% sub_combined$PIT
    if(fishID %in% sub_combined$PIT){
      visit_count<-visit_count+1
    }
    }
    superfish$nvisit[i]<-visit_count
}
  return(superfish)
}

Results

Bipartite networks for Lake Parley 2020 season

#Bipartite network for boat ramp site
edgelist<-data.frame(PIT=merged_all_sites$PIT, Reader=merged_all_sites$Site)

A <- spMatrix(nrow=length(unique(edgelist$PIT)),
        ncol=length(unique(edgelist$Reader)),
        i = as.numeric(factor(edgelist$PIT)),
        j = as.numeric(factor(edgelist$Reader)),
        x = rep(1, length(as.numeric(edgelist$PIT))) )
row.names(A) <- levels(factor(edgelist$PIT))
colnames(A) <- levels(factor(edgelist$Reader))



#using carp data
bi<-graph.incidence(A, mode="all") #undirected, named graph that is bipartite
pr<-bipartite.projection(bi, multiplicity=TRUE) 
#co-membership network of nodes ($proj1), or a network of groups that share members ($proj2)

# get.adjacency(pr$proj1,sparse=FALSE,attr="weight")
h<-pr$proj2
l3 <-layout.circle(h)
plot.igraph(h,vertex.label=V(h)$name,layout=l3,edge.width=E(h)$weight/10^8, edge.color="black")

# visweb(t(A), circles=TRUE,  boxes=FALSE,  labsize=1, circle.max=1.5, text="no", plotsize= 18)
plotweb(as.matrix(A), text.rot=90, labsize=1, low.lablength = 6, text.low.col="black")

Number of visits through time

Below we show the number of visits through time for eight sites in Lake Parley for the 2020 season from June 23- September 9, 2020. The two blue lines in the histograms below indicate the dates of baiting (July 3, 2020) and adding nets plus the first capture event respectively (July 23, 2020).

#histogram of number of unique visitors by sampling day
start_date<-as.POSIXct("2020-06-23 12:00:00 CDT")
end_date<-as.POSIXct("2020-09-09 23:59:59 CDT")

visits_through_time<-function(start_date,end_date, merged_df, site_name){

num_days<-as.numeric(end_date-start_date)
sampling_df <-data.frame(day=1:num_days, uniqueID=rep(NA, num_days))
degree_list<-NULL

for(i in 1:num_days){
  sub_combined<-merged_df %>% filter(DATETIME >= as.POSIXct(start_date+ (i-1)*86400, tz="CDT") & DATETIME <= as.POSIXct(start_date + i*86400, tz="CDT")) #86400 seconds/24 hours
  # print(sub_combined)
  sampling_df$uniqueID[i]<-length(unique(sub_combined$PIT))
  PITs<-unique(sub_combined$PIT)
  if(length(PITs)>0){
  carp.degree<-data.frame(PITs, numsites= rep(NA, length(PITs)))
      for(j in 1:length(PITs)){
        carp.degree$numsites[j]<-length(unique(sub_combined$Antenna[which(sub_combined$PIT==PITs[j])]))
      }
  degree_list[[i]]<-carp.degree
  }
}

sampling_df$uniqueID<-as.numeric(sampling_df$uniqueID)
barplot(uniqueID ~day, data=sampling_df, xlab="Sampling Day", ylab="Number of unique carp", main=paste(site_name, expression("Number of unique carp detected per day")))
abline (v = 10, col="blue") #corn/bait, July 3, 2020
abline (v= 30, col="blue") #net deployed and capture attempt, July 23, 2020
}

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s1, site_name="Site 1:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s2, site_name="Site 2:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s3, site_name="Site 3:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s4, site_name="Site 4:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s5, site_name="Site 5:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s6, site_name="Site 6:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s7, site_name="Site 7:")

visits_through_time(start_date=start_date, end_date=end_date, merged_df=merged_s8, site_name="Site 8:")

Visits during baiting period

In these plots, individual visits of unique carp are shown through time at each site.

#subdivide data based on phases of experiment
#corn/bait, July 3, 2020
#first net deployment and capture attempt, July 23, 2020- 10 AM

# min(merged_br$DATETIME)
# max(merged_br$DATETIME)

s1_corn<- merged_s1 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s2_corn<- merged_s2 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s3_corn<- merged_s3 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s4_corn<- merged_s4 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s5_corn<- merged_s5 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s6_corn<- merged_s6 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s7_corn<- merged_s7 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
s8_corn<- merged_s8 %>% filter(DATETIME >=as.POSIXct("2020-07-03 12:00:00", tz="CDT") & DATETIME < as.POSIXct("2020-07-23 12:00:00", tz="CDT"))
#time course of individual carp visits to each site during baiting period (July 3rd- July 23rd 2020)
ggplot(data=s1_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 1")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s2_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 2")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s3_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 3")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s4_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 4")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s5_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 5")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s6_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 6")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s7_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 7")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

ggplot(data=s8_corn, aes(DATETIME, as.factor(PIT))) +
  geom_point(alpha=0.5)+
  scale_x_datetime(date_breaks="1 day")+
  ggtitle("Site 8")+
  xlab("Date")+
  ylab("Carp PIT ID")+
  theme(axis.text.x=element_text(angle=90), axis.text.y=element_text(size=6), axis.title = element_text(size=10))

Superfeeders and feeding preference

Below we show histograms for the number of sampling days for all eight sites.

start_date<-as.POSIXct("2020-07-03 12:00:00 CDT") #rounded to nearest midday
end_date<-as.POSIXct("2020-07-23 12:00:00 CDT")

num_days<-as.numeric(end_date-start_date)

#Histogram of the number of sampling days each fish visits Site 1
superfish_s1<-calc_nvisit(df=s1_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s1,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 1: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 2
superfish_s2<-calc_nvisit(df=s2_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s2,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 2: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 3
superfish_s3<-calc_nvisit(df=s3_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s3,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 3: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 4
superfish_s4<-calc_nvisit(df=s4_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s4,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 4: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 5
superfish_s5<-calc_nvisit(df=s5_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s5,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 5: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 6
superfish_s6<-calc_nvisit(df=s6_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s6,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 6: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 7
superfish_s7<-calc_nvisit(df=s7_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s7,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 7: Baiting Period")

#Histogram of the number of sampling days each fish visits Site 8
superfish_s8<-calc_nvisit(df=s8_corn, num_days=num_days, start_date=start_date)
ggplot(superfish_s8,aes(nvisit)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="light blue") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of sampling days")+
  ylab("Number of unique fish")+
  ggtitle("Site 8: Baiting Period")

We also combine visits by site and show the relative proportion of visits to each site per fish. We highlight “superfeeders”–those fish accounting for the top 20% of total visits to the left of the vertical black line.

#Combined plot showing site preference and number of detected visits by fish
superfish_s1$Site<-"Site1"
superfish_s2$Site<-"Site2"
superfish_s3$Site<-"Site3"
superfish_s4$Site<-"Site4"
superfish_s5$Site<-"Site5"
superfish_s6$Site<-"Site6"
superfish_s7$Site<-"Site7"
superfish_s8$Site<-"Site8"

all_superfish <- rbind(superfish_s1,superfish_s2,superfish_s3,superfish_s4,superfish_s5,superfish_s6,superfish_s7,superfish_s8, Fill=NA)
all_superfish<-na.omit(all_superfish)
superfish_wide <- spread(all_superfish, Site, nvisit) #convert from long to wide data format
# superfish_wide <- subset(superfish_wide,select = fishID:Site8)
superfish_wide[is.na(superfish_wide)] <- 0 #convert NAs to zero
superfish_wide$total <- rowSums(superfish_wide[2:9])
superfish_wide<-superfish_wide[order(superfish_wide$total, decreasing=TRUE),]

all_superfish$fishID <- factor(all_superfish$fishID, levels = superfish_wide$fishID[order(superfish_wide$total, decreasing=TRUE)])

top20<-round(length(levels(all_superfish$fishID))*.2)
cutoff<-levels(all_superfish$fishID)[top20]

ggplot(all_superfish, aes(fill=Site, y=nvisit, x=fishID))+ 
    geom_col()+
    theme(axis.text.x = element_text(size=6, angle = 90, hjust = 1))+
    scale_fill_viridis(discrete=TRUE, name = "Site", labels = c("Site 1", "Site 2", "Site 3", "Site 4", "Site 5", "Site 6", "Site 7", "Site 8"))+
    labs(y= "Visits", x = "Carp ID")+
    geom_vline(xintercept=which(levels(all_superfish$fishID) == cutoff))

Average group size and feeding bout duration

Below we show histograms for average group size during feeding for sites 1-8, as well as all sites combined, as detected by Gaussian mixed models (GMM).

events_s1 <- readRDS("~/Carp/2020/events_s1.RDS")
events_s2 <- readRDS("~/Carp/2020/events_s2.RDS")
events_s3 <- readRDS("~/Carp/2020/events_s3.RDS")
events_s4 <- readRDS("~/Carp/2020/events_s4.RDS")
events_s5 <- readRDS("~/Carp/2020/events_s5.RDS")
events_s6 <- readRDS("~/Carp/2020/events_s6.RDS")
events_s7 <- readRDS("~/Carp/2020/events_s7.RDS")
events_s8 <- readRDS("~/Carp/2020/events_s8.RDS")
events_parley2020 <- readRDS("~/Carp/2020/events_parley2020.RDS")

gbi_s1 <- readRDS("~/Carp/2020/gbi_s1.RDS")
gbi_s2 <- readRDS("~/Carp/2020/gbi_s2.RDS")
gbi_s3 <- readRDS("~/Carp/2020/gbi_s3.RDS")
gbi_s4 <- readRDS("~/Carp/2020/gbi_s4.RDS")
gbi_s5 <- readRDS("~/Carp/2020/gbi_s5.RDS")
gbi_s6 <- readRDS("~/Carp/2020/gbi_s6.RDS")
gbi_s7 <- readRDS("~/Carp/2020/gbi_s7.RDS")
gbi_s8 <- readRDS("~/Carp/2020/gbi_s8.RDS")
gbi_parley2020 <- readRDS("~/Carp/2020/gbi_parley2020.RDS")
#histogram of group size for Site 1
s1_gs<-data.frame(bout=1:nrow(gbi_s1), groupsize=rowSums(gbi_s1))
s1_gs$dur<-events_s1$End-events_s1$Start
s1_gs$dur_min<-s1_gs$dur/60
head(s1_gs)
##   bout groupsize  dur   dur_min
## 1    1         2  346  5.766667
## 2    2         2 1211 20.183333
## 3    3         2   90  1.500000
## 4    4         2  486  8.100000
## 5    5         2  378  6.300000
## 6    6         5 2662 44.366667
ggplot(s1_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 1: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s1_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 1: Duration of feeding bout (min)")

mean(s1_gs$dur_min)
## [1] 5.977423
median(s1_gs$dur_min)
## [1] 4.1
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s1_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 1")

#histogram of group size for Site 2
s2_gs<-data.frame(bout=1:nrow(gbi_s2), groupsize=rowSums(gbi_s2))
s2_gs$dur<-events_s2$End-events_s2$Start
s2_gs$dur_min<-s2_gs$dur/60
head(s2_gs)
##   bout groupsize  dur   dur_min
## 1    1         2  703 11.716667
## 2    2         2 1943 32.383333
## 3    3         2  107  1.783333
## 4    4         2  328  5.466667
## 5    5         2  528  8.800000
## 6    6         2  554  9.233333
ggplot(s2_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 2: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s2_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 2: Duration of feeding bout (min)")

mean(s2_gs$dur_min)
## [1] 4.597763
median(s2_gs$dur_min)
## [1] 2.966667
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s2_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 2")

#histogram of group size for Site 3
s3_gs<-data.frame(bout=1:nrow(gbi_s3), groupsize=rowSums(gbi_s3))
s3_gs$dur<-events_s3$End-events_s3$Start
s3_gs$dur_min<-s3_gs$dur/60
head(s3_gs)
##   bout groupsize dur   dur_min
## 1    1         2 302 5.0333333
## 2    2         2  41 0.6833333
## 3    3         4 578 9.6333333
## 4    4         2 289 4.8166667
## 5    5         3 487 8.1166667
## 6    6         2 294 4.9000000
ggplot(s3_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 3: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s3_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 3: Duration of feeding bout (min)")

mean(s3_gs$dur_min)
## [1] 4.761374
median(s3_gs$dur_min)
## [1] 3.566667
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s3_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 3")

#histogram of group size for Site 4
s4_gs<-data.frame(bout=1:nrow(gbi_s4), groupsize=rowSums(gbi_s4))
s4_gs$dur<-events_s4$End-events_s4$Start
s4_gs$dur_min<-s4_gs$dur/60
head(s4_gs)
##   bout groupsize dur  dur_min
## 1    1         2 305 5.083333
## 2    2         2 524 8.733333
## 3    3         2 223 3.716667
## 4    4         2 304 5.066667
## 5    5         2 390 6.500000
## 6    6         2 250 4.166667
ggplot(s4_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 4: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s4_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 4: Duration of feeding bout (min)")

mean(s4_gs$dur_min)
## [1] 5.592101
median(s4_gs$dur_min)
## [1] 4.8
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s4_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 4")

#histogram of group size for Site 5
s5_gs<-data.frame(bout=1:nrow(gbi_s5), groupsize=rowSums(gbi_s5))
s5_gs$dur<-events_s5$End-events_s5$Start
s5_gs$dur_min<-s5_gs$dur/60
head(s5_gs)
##   bout groupsize  dur   dur_min
## 1    1         2 1300 21.666667
## 2    2         2  326  5.433333
## 3    3         2  116  1.933333
## 4    4         2  227  3.783333
## 5    5         3  758 12.633333
## 6    6         2 2786 46.433333
ggplot(s5_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 5: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s5_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 5: Duration of feeding bout (min)")

mean(s5_gs$dur_min)
## [1] 8.76577
median(s5_gs$dur_min)
## [1] 4.15
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s5_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 5")

#histogram of group size for Site 6
s6_gs<-data.frame(bout=1:nrow(gbi_s6), groupsize=rowSums(gbi_s6))
s6_gs$dur<-events_s6$End-events_s6$Start
s6_gs$dur_min<-s6_gs$dur/60
head(s6_gs)
##   bout groupsize dur   dur_min
## 1    1         2  41 0.6833333
## 2    2         6 345 5.7500000
## 3    3         4 440 7.3333333
## 4    4         3 168 2.8000000
## 5    5         3 570 9.5000000
## 6    6         3 420 7.0000000
ggplot(s6_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 6: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s6_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 6: Duration of feeding bout (min)")

mean(s6_gs$dur_min)
## [1] 4.195252
median(s6_gs$dur_min)
## [1] 3.033333
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s6_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 6")

#histogram of group size for Site 7
s7_gs<-data.frame(bout=1:nrow(gbi_s7), groupsize=rowSums(gbi_s7))
s7_gs$dur<-events_s7$End-events_s7$Start
s7_gs$dur_min<-s7_gs$dur/60
head(s7_gs)
##   bout groupsize dur   dur_min
## 1    1         2 178 2.9666667
## 2    2         2 196 3.2666667
## 3    3         2  59 0.9833333
## 4    4         3 312 5.2000000
## 5    5         3 418 6.9666667
## 6    6         3 248 4.1333333
ggplot(s7_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 7: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s7_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 7: Duration of feeding bout (min)")

mean(s7_gs$dur_min)
## [1] 3.48897
median(s7_gs$dur_min)
## [1] 1.983333
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s7_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 7")

#histogram of group size for Site 8
s8_gs<-data.frame(bout=1:nrow(gbi_s8), groupsize=rowSums(gbi_s8))
s8_gs$dur<-events_s8$End-events_s8$Start
s8_gs$dur_min<-s8_gs$dur/60
head(s8_gs)
##   bout groupsize dur  dur_min
## 1    1         2 198 3.300000
## 2    2         2 421 7.016667
## 3    3         2 162 2.700000
## 4    4         3 430 7.166667
## 5    5         2 140 2.333333
## 6    6         2 240 4.000000
ggplot(s8_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Site 8: Average group size")

#histogram of feeding bout duration for Site 1
ggplot(s8_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Site 8: Duration of feeding bout (min)")

mean(s8_gs$dur_min)
## [1] 4.508432
median(s8_gs$dur_min)
## [1] 3.583333
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(s8_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Site 8")

#histogram of group size for all of Lake Parley
parley2020_gs<-data.frame(bout=1:nrow(gbi_parley2020), groupsize=rowSums(gbi_parley2020))
parley2020_gs$dur<-events_parley2020$End-events_parley2020$Start
parley2020_gs$dur_min<-parley2020_gs$dur/60
head(parley2020_gs)
##   bout groupsize  dur   dur_min
## 1    1         2 1139 18.983333
## 2    2         3  443  7.383333
## 3    3         2  270  4.500000
## 4    4         2  751 12.516667
## 5    5         2  174  2.900000
## 6    6         2  494  8.233333
ggplot(parley2020_gs, aes(groupsize)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous(breaks=0:20)+
  xlab("Number of carp per feeding bout")+
  ylab("Number of times observed")+
  ggtitle("Lake Parley: Average group size")

#histogram of feeding bout duration for all of Lake Parley
ggplot(parley2020_gs, aes(dur_min)) + 
geom_histogram(binwidth=1,boundary=-0.5, color="black", fill="darkcyan") + 
scale_x_continuous()+
  xlab("Duration of feeding bout (min)")+
  ylab("Number of times observed")+
  ggtitle("Lake Parley: Duration of feeding bout (min)")

mean(parley2020_gs$dur_min)
## [1] 5.111328
median(parley2020_gs$dur_min)
## [1] 3.066667
#QQ plot of group size vs. bout duration Boat Ramp (separate antenna)
ggplot(parley2020_gs, aes(x=dur_min, y=groupsize)) +
  geom_point(size=2, shape=19, alpha=0.6)+
  xlab("Bout duration (min)")+
  ylab("Group size")+
  ggtitle("Lake Parley")

GMM Networks

Here we show plots of individual networks derived from sites 1-8 and Lake Parley as a whole. This code also produced histograms of weighted degree which reflects both the number of edges (contacts) of a given node, but also the frequency or intensity of those interactions.

Site 1 Networks

library(asnipe)
library(igraph)

# Network of Site 1
identities<-colnames(gbi_s1)
network <- matrix(0, ncol(gbi_s1), ncol(gbi_s1))

network<- get_network(gbi_s1, data_format="GBI",
association_index="HWI", identities=identities, times=events_s1$midpoint, start_time=min(events_s1$midpoint), end_time=max(events_s1$midpoint))
## Generating  33  x  33  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 1 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 1")

network_stats_s1<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s1$node_names<-rownames(network_stats_s1)

s1_long<-gather(network_stats_s1, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s1, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 1: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 2 Networks

library(asnipe)
library(igraph)

# Network of Site 2
identities<-colnames(gbi_s2)
network <- matrix(0, ncol(gbi_s2), ncol(gbi_s2))

network<- get_network(gbi_s2, data_format="GBI",
association_index="HWI", identities=identities, times=events_s2$midpoint, start_time=min(events_s2$midpoint), end_time=max(events_s2$midpoint))
## Generating  50  x  50  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 2 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 2")

network_stats_s2<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s2$node_names<-rownames(network_stats_s2)

s2_long<-gather(network_stats_s2, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s2, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 2: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 3 Networks

library(asnipe)
library(igraph)

# Network of Site 3
identities<-colnames(gbi_s3)
network <- matrix(0, ncol(gbi_s3), ncol(gbi_s3))

network<- get_network(gbi_s3, data_format="GBI",
association_index="HWI", identities=identities, times=events_s3$midpoint, start_time=min(events_s3$midpoint), end_time=max(events_s3$midpoint))
## Generating  41  x  41  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 3 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 3")

network_stats_s3<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s3$node_names<-rownames(network_stats_s3)

s3_long<-gather(network_stats_s3, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s3, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 3: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 4 Networks

library(asnipe)
library(igraph)

# Network of Site 4
identities<-colnames(gbi_s4)
network <- matrix(0, ncol(gbi_s4), ncol(gbi_s4))

network<- get_network(gbi_s4, data_format="GBI",
association_index="HWI", identities=identities, times=events_s4$midpoint, start_time=min(events_s4$midpoint), end_time=max(events_s4$midpoint))
## Generating  30  x  30  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 4 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 4")

network_stats_s4<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s4$node_names<-rownames(network_stats_s4)

s4_long<-gather(network_stats_s4, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s4, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 4: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 5 Networks

library(asnipe)
library(igraph)

# Network of Site 5
identities<-colnames(gbi_s5)
network <- matrix(0, ncol(gbi_s5), ncol(gbi_s5))

network<- get_network(gbi_s5, data_format="GBI",
association_index="HWI", identities=identities, times=events_s5$midpoint, start_time=min(events_s5$midpoint), end_time=max(events_s5$midpoint))
## Generating  34  x  34  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 5 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 5")

network_stats_s5<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s5$node_names<-rownames(network_stats_s5)

s5_long<-gather(network_stats_s5, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s5, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 5: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 6 Networks

library(asnipe)
library(igraph)

# Network of Site 6
identities<-colnames(gbi_s6)
network <- matrix(0, ncol(gbi_s6), ncol(gbi_s6))

network<- get_network(gbi_s6, data_format="GBI",
association_index="HWI", identities=identities, times=events_s6$midpoint, start_time=min(events_s6$midpoint), end_time=max(events_s6$midpoint))
## Generating  42  x  42  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 6 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 6")

network_stats_s6<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s6$node_names<-rownames(network_stats_s6)

s6_long<-gather(network_stats_s6, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s6, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 6: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 7 Networks

library(asnipe)
library(igraph)

# Network of Site 7
identities<-colnames(gbi_s7)
network <- matrix(0, ncol(gbi_s7), ncol(gbi_s7))

network<- get_network(gbi_s7, data_format="GBI",
association_index="HWI", identities=identities, times=events_s7$midpoint, start_time=min(events_s7$midpoint), end_time=max(events_s7$midpoint))
## Generating  37  x  37  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 7 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 7")

network_stats_s7<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s7$node_names<-rownames(network_stats_s7)

s7_long<-gather(network_stats_s7, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s7, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 7: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Site 8 Networks

library(asnipe)
library(igraph)

# Network of Site 8
identities<-colnames(gbi_s8)
network <- matrix(0, ncol(gbi_s8), ncol(gbi_s8))

network<- get_network(gbi_s8, data_format="GBI",
association_index="HWI", identities=identities, times=events_s8$midpoint, start_time=min(events_s8$midpoint), end_time=max(events_s8$midpoint))
## Generating  37  x  37  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Site 8 (SRI)")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label=NA, vertex.label.color="black", main="Site 8")

network_stats_s8<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_s8$node_names<-rownames(network_stats_s8)

s8_long<-gather(network_stats_s8, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_s8, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Site 8: Degree and betweenness")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Lake Parley 2020 Networks

# Network of Parley 
identities<-colnames(gbi_parley2020)
network <- matrix(0, ncol(gbi_parley2020), ncol(gbi_parley2020))

network<- get_network(gbi_parley2020, data_format="GBI",
association_index="HWI", identities=identities, times=events_parley2020$midpoint, start_time=min(events_parley2020$midpoint),     end_time=max(events_parley2020$midpoint))
## Generating  133  x  133  matrix
rownames(network)<-identities
colnames(network)<-identities

net <- graph.adjacency(network, mode="undirected", diag=FALSE, weighted=TRUE)
deg_weighted <- graph.strength(net)
hist(deg_weighted, main="Histogram of Weighted Degree for Parley Network")

V(net)$size <- deg_weighted*10
plot(net, vertex.color="cyan4", vertex.label= NA, main="Parley: Combined Antenna")

network_stats_parley2020<-data.frame(nodes=as.character(V(net)), weighted_degree=igraph::graph.strength(net), degree=igraph::degree(net), betweenness=igraph::betweenness(net),transitivity=
igraph::transitivity(net, type="localundirected"))
network_stats_parley2020$node_names<-rownames(network_stats_parley2020)

parley2020_long<-gather(network_stats_parley2020, key=network_metric, value=value,
weighted_degree:transitivity)

ggplot(network_stats_parley2020, aes(x=weighted_degree, y=betweenness))+
  geom_point(aes(size=transitivity, alpha=0.5))+
  xlab("Weighted degree")+
  ylab("Betweenness")+
  ggtitle("Parley (all sites)")+
  geom_text(aes(label=node_names),hjust=0, vjust=0)

Boxplots of individual network measures

Here we plot boxplots of individual unweighted degree, overall association strength (weighted degree) and betweenness for networks of Sites 1-8 and Lake Parley as a whole.

# boxplots of individual network measures 
ggplot(s1_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 1")+ #
  xlab("Network metric")

ggplot(s2_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 2")+ #
  xlab("Network metric")

ggplot(s3_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 3")+ #
  xlab("Network metric")

ggplot(s4_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 4")+ #
  xlab("Network metric")

ggplot(s5_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 5")+ #
  xlab("Network metric")

ggplot(s6_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 6")+ #
  xlab("Network metric")

ggplot(s7_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 7")+ #
  xlab("Network metric")

ggplot(s8_long, aes(x=network_metric, y=value)) + 
  geom_boxplot()+ 
  ggtitle("Network metrics for Site 8")+ #
  xlab("Network metric")

ggplot(parley2020_long, aes(x=network_metric, y=value)) + #
  geom_boxplot()+ 
  ggtitle("Network metrics for Lake Parley")+ #
  xlab("Network metric") 

s1_long$site<-rep("Site 1", nrow(s1_long))
s2_long$site<-rep("Site 2", nrow(s2_long))
s3_long$site<-rep("Site 3", nrow(s3_long))
s4_long$site<-rep("Site 4", nrow(s4_long))
s5_long$site<-rep("Site 5", nrow(s5_long))
s6_long$site<-rep("Site 6", nrow(s6_long))
s7_long$site<-rep("Site 7", nrow(s7_long))
s8_long$site<-rep("Site 8", nrow(s8_long))

parley2020_long$site<-rep("Parley", nrow(parley2020_long))

all_sites_long<-rbind(s1_long, s2_long, s3_long, s4_long, s5_long, s6_long, s7_long, s8_long, parley2020_long)

ggplot(all_sites_long, aes(x=site, y=value)) + #
  geom_boxplot()+ 
  facet_grid(network_metric~., scales = "free")+
  ggtitle("Comparing network metrics across sites")+ #
  xlab("Site") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Reproducibility

Sys.time()
## [1] "2021-03-18 10:59:59 EDT"
git2r::repository()
## Local:    master /research-home/lwhite/Carp
## Remote:   master @ origin (https://github.com/whit1951/Carp.git)
## Head:     [7d9cb9d] 2020-07-01: Current html output
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.7 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] asnipe_1.1.12        viridis_0.5.1        viridisLite_0.3.0   
##  [4] tidyr_1.1.3          gganimate_1.0.5      ggplot2_3.3.3       
##  [7] data.table_1.14.0    bipartite_2.15       sna_2.6             
## [10] network_1.16.1       statnet.common_4.4.1 vegan_2.5-6         
## [13] lattice_0.20-41      permute_0.9-5        igraph_1.2.6        
## [16] Matrix_1.3-2         lubridate_1.7.10     stringr_1.4.0       
## [19] dplyr_1.0.5         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        prettyunits_1.1.1 assertthat_0.2.1  digest_0.6.27    
##  [5] utf8_1.2.1        rle_0.9.2         R6_2.5.0          evaluate_0.14    
##  [9] coda_0.19-4       spam_2.5-1        highr_0.8         pillar_1.5.1     
## [13] rlang_0.4.10      progress_1.2.2    jquerylib_0.1.3   rmarkdown_2.7    
## [17] labeling_0.4.2    splines_4.0.3     munsell_0.5.0     compiler_4.0.3   
## [21] xfun_0.22         pkgconfig_2.0.3   mgcv_1.8-34       htmltools_0.5.1.1
## [25] tidyselect_1.1.0  gridExtra_2.3     tibble_3.1.0      codetools_0.2-18 
## [29] fansi_0.4.2       crayon_1.4.1      withr_2.4.1       MASS_7.3-53.1    
## [33] grid_4.0.3        nlme_3.1-152      jsonlite_1.7.2    gtable_0.3.0     
## [37] lifecycle_1.0.0   DBI_1.1.1         git2r_0.28.0      magrittr_2.0.1   
## [41] scales_1.1.1      stringi_1.5.3     debugme_1.1.0     farver_2.1.0     
## [45] bslib_0.2.4       ellipsis_0.3.1    generics_0.1.0    vctrs_0.3.6      
## [49] tools_4.0.3       glue_1.4.2        tweenr_1.0.1      purrr_0.3.4      
## [53] maps_3.3.0        hms_1.0.0         fields_10.3       parallel_4.0.3   
## [57] yaml_2.2.1        colorspace_2.0-0  cluster_2.1.1     dotCall64_1.0-0  
## [61] knitr_1.31        sass_0.3.1

References

Damien R. Farine (2019). asnipe: Animal Social Network Inference and Permutations for Ecologists. R package version 1.1.12. https://CRAN.R-project.org/package=asnipe